Phase field crystal modeling of early stage precipitation and clustering 
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A phase field crystal model is used to investigate the mechanisms of formation and growth of 
early clusters in quenched/aged dilute binary alloys, a phenomenon typically outside the scope of 
molecular dynamics time scales. We show that formation of early sub-critical clusters is triggered 
by the stress relaxation effect of quenched-in defects, such as dislocations, on the energy barrier and 
the critical size for nucleation. In particular, through analysis of system energetics, we demonstrate 
that the growth of sub-critical clusters into overcritical sizes occurs due to the fact that highly 
strained areas in the lattice locally reduce or even eliminate the free energy barrier for a first-order 
transition. 



I. INTRODUCTION 

Clustering in metal alloys is known as the very early 
stage of first-order transformations within a bulk crys- 
tal, which largely influences the mechanical properties 
of quenched/aged materials. Coherent with the matrix, 
small clusters of solute atoms have a signiflcantly lower 
nucleation barrier than the terminal second-phase of a 
completely different crystal structure. Early clusters nu- 
cleate and grow in size forming the so-called CP-zones 
known to be a metastable precursor of the equilibrium 
phasd^. The formation and growth mechanisms of the 
early clusters are poorly understood presently due to 
the lack of direct atomistic observations of structural 
changes during the transition process. However, inspired 
by the observations made on the quenched structures us- 
ing tr ansm ission electron microscopy (TEMP^, 

3D atom 

prob^^lsMH^ and positron annihilatiorPSMl techniques, 
the formation of early clusters has been empirically asso- 
ciated with so-called quenched-in defects. Formed within 
the bulk crystal upon quenching, excess vacancies and/or 
dislocations loops have been presumed to decrease the 
ene rgy barr ier for nucleation facilitating cluster forma- 
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Understanding solute clustering mechanisms is of cru- 
cial importance to design an effective age-hardening pro- 
cess producing desired mechanical properties in alloys^. 
To our knowledge, no systematic study of the clustering 
mechanisms has been done using atomic-scale simulation 
methods such as molecular dynamics (MD) and Monte 
Carlo (MC) simulations, due to their main restrictions of 
accessing the relevant time scales of diffusional transfor- 
mations. Dynamical calculations using classical density 
functional theory (CDFT) are also inefficient due to the 
high spatial resolution required to resolve the sharp den- 
sity spikes in solid phases^. 

The recently developed phase field crystal (PFC) 
methocP^HUl has shown promise for simulating structural 



transformations on diffusive time scales. This new for- 
malism carries the essential physics of CDFT without 
the need to resolve the sharp atomic density peaks. In 
the most recent PFC formalism developed by Greenwood 
et aL^i^ltl^, various crystal symmetries can be easily stabi- 
lized by construction of relevant correlation kernels. This 
approach preserves the numerical efficiency of the original 
PFC model and is able to dynamically simulate the pre- 
cipitation of solid phases within a parent phase of differ- 
ent crystal symmetry^^ and/or chemical compositiorP^. 

This letter proposes a new approach to study the clus- 
tering phenomenon that relies on atomic-scale simula- 
tions using the previously developed alloy PFC model of 
ref.13. We explore the formation and growth mechanisms 
of early clusters in a quenched bulk lattice of a supersatu- 
rated Al-Cu alloy initially containing quenched-in defects 
such as dislocations. 

II. MODEL STRUCTURE 

We start with the free energy functional in the binary 
PFC modeli^, 

(n -I- l)Af,„„ -^nj dr'C:ff{\r - r'\)n' + alVcpjdr 

(1) 

where n and c represent reduced dimensionless atomic 
number density and solute concentration fields, respec- 
tively, rj and x are coefficients added to fit the ideal 
energy to a polynomial expansion (?7 = X = 1 describes 
a Taylor series expansion of the bulk free energy around 
the reference density) and 

Af™„ =w{cln(-) + (l-c)ln(-^^)} (2) 
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represents the energy density associated with the entropy 
of mixing. The coefhcient uj is introduced to fit the en- 
tropic energy away from the reference composition cq. 
The parameter a is a coefficient (taken as 1 in this study). 
These parameters are discussed further in ref.^'*'. 

For a binary alloy, Greenwood et al.^^ introduced the 
correlation function 



c:ff^x^{c)c^^ + X2{c)a 
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, where Xi{c) = l-3c^ + 2c^ andX2(c) = l-3(l-c)2 + 
2(1 — c)^. and are correlation functions rep- 

resenting, respectively, contributions to the excess free 
energy for the situations where A atoms are in the pre- 
ferred crystalline network of B atoms and B atoms which 
are in a structure preferred by A atoms. The correla- 
tion functions Cf{k) are defined to have reciprocal space 
peaks (i.e. fcj, corresponding to the inverse of interpla- 
nar spacings) determined by the main families of planes 
in the equilibrium crystal unit cell structure for the i*'' 
component. Each peak is represented by the following 
Gaussian form of width a^, modulated for temperature 
by a Debye- Waller prefactor which accounts for an effec- 
tive transition temperature (JmJ^- 
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The equations of motion of the total density and con- 
centration fields follow dissipative dynamicJ^^. The total 
mass density and total reference density per unit volume 
are defined as p = pA + Pb and p° = P°a + P°b, respec- 
tively. Thus, the equations of motion can be written for 
n(= p/ p° — 1) and c(= 



Pb/p) as ^ 



V.{M„V(^)}- 



?7„(a,i) and ff = V.{MeV(^)} + 77,((t, i), respec- 
tiveljEI. Mn and Mc are dimensionless kinetic mobility 
parameters (equal to 1 in this study). ?7n(o', t) and ?7c(ct, t) 
are stochastic noise variables subsuming the role of fast 
atomic vibrations in density and concentration fields, re- 
spectively. 



III. RESULTS 

A. Phase diagram reconstruction 

To examine the equilibrium properties of this binary 
PFG model for a 2D Al-Cu system, we construct the 
phase diagram for the coexistence of two square phases. 
The coexistence lines between the respective phases are 
obtained by a common tangent construction of the free 
energy curves of solid and liquid at the reference density 
[n = 0). Following Greenwood et al.'i^, the free energy 
curves of the square phases are calculated using the two- 
mode approximation of the density fields which is defined 

by 
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, where the subscript i denotes a particular solid phase 
with a lattice spacing a^, and the index j counts the Ni 
modes of the i-phase. Aj is the amplitude of mode j and 
I is the index over the group of reciprocal space peaks 
corresponding to mode j, Nj. Accordingly, fcjj is the 
reciprocal lattice vector normalized to a lattice spacing 
of 1, corresponding to each index I in the family j. The 
free energy curve for each phase can be calculated as a 
function of the composition c by substituting the above 
density field approximation into Eq. ([T]) and integrating 
over the unit cell. The resulting crystal free energy is 
then minimized for the amplitudes Aj. For the liquid 
phase, the amplitude Aj is set to zero and the density 
is considered as the reference density {n — {)). A more 
detailed description of this methodology is provided in 
ref.T^. 
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FIG. 1. The constructed phase diagram for a square-square 
system with the inset showing the Al-rich side of the exper- 
imental phase diagram of Al-Cu system taken from Ref.'^. 
The parameters for ideal free energy contribution were rj = 1.4 
and X = 1) while uj — 0.005 and Co = 0.5 for entropy of 
mixing. Widths of the correlations peaks are Oua; = 2.4, 
ceioAi = \Ploi\\Ai (the required ratio to introduce isotropic 
elastic constants in an square phasePl), 

One = 2.4 and 

aioe ~ V^aiio- The peak positions for pure Al correspond 
to kiiAi = 2tv, kioAi = V^kiiAi, kiw = (81/38)7r and 
kioe = \/2fciie- The effective transition temperatures are 
set to OMwAi — 0.55, omwai — 0.55, omiib = 0.55 and 
o"A/ioe = 0.55; The concentration c is rescaled considering the 
Cu-content in the S-phase. 

In the Al-rich side of the experimental Al-Cu phase 
diagram, shown in the inset of Fig. [l] there is a eutectic 
transition between the Al-rich a-fcc phase and an in- 
termediate phase 9 (containing « 32.5at.% Cu) with a 
tetragonal crystal structure. For 2D simulations, in order 
to approximate these equilibrium properties, we recon- 
struct the binary phase diagram of Al and 0, both with 
a square symmetry but differing in Cu-content. The lat- 
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tice constant (and thus the reciprocal space peaks) of is 
approximated by interpolating between those of Pure Al 
and Cu. The solid phase free energy is calculated with 
a variable lattice constant weighted by concentration c 
using the interpolation functions Xi and X2. The poly- 
nomial fitting parameters in Eq. ([T]) (namely 77, x a-nd uj) 
and width of various peaks {oij) in the correlation kernel 
C^j are then chosen so as to obtain the same composi- 
tions for a-phase solubility limit and eutectic point as 
those in the experimental phase diagram. 

B. Simulation of clustering 

With the equilibrium properties obtained above, sim- 
ulations of clustering were performed on a rectangu- 
lar mesh with grid spacing dx = 0.125 and time step 
dt = \. Considering the lattice parameter of 1, each 
atomic spacing is resolved by 8 mesh spacings. The dy- 
namical equations were solved semi-implicitly in Fourier 
space for higher efficiency. The initial conditions were 
chosen to study the proposed dominant role of quenched- 
in dislocation-type defects in the bulk crystal during the 
early stage precipitation in dilut e Al-Cu a lloys quenched 
from a solutionizing temperatur c!^ l ^ l ^ l ^° l ^° l According to 
this hypothesis, dislocation loops, generated by excess va- 
cancies, are responsible for local lattice distortions facil- 
itating segregation and diffusion of Cu-atoms, while also 
driving the s yste m towards a more thermodynamically- 
stable state^iEIini Therefore, as initial conditions, we use 
a crystal lattice of uniform composition distorted by in- 
troducing dislocations. 

PEC simulation is performed for quench/aging of Al- 
2at.%Cu from the solutionizing temperature of ct = 0.17 
to fj = 0.04 with the initial conditions shown in Fig.[2][b). 
During the simulation, first, small clusters form with a 
slightly higher Cu-content than that of the matrix. As 
time progresses, some of these clusters shrink in size and 
concentration and a few get stabilized (e.g. the cluster 
shown in Eig. [2|^a)). In contrast, as expected, quenching 
the same initial structure from the solutionizing tempera- 
ture of (T = 0.17 to a temperature within the single-phase 
Sq-Al region, i.e., a ~ 0.16, leads to complete removal of 
distortion. 



IV. DISCUSSION 

A. Evolution of clusters 

The dislocation-induced cluster structure shown in 
Fig. [2]ja) is consistent with TEM observations in 
Al-1.7at.%CiPl and Al-l.lat.%Cu-0.5at.%Mg alloys^^ 
where dislocation loops appear in the bulk lattice of the 
quenched structures. Using resistometric measurements 
and T EM t echniques for Al-1.2at.%Si alloy, Ozawa and 
KimurgPlil have associated the formation of dislocation 
(or vacancy) loops upon quenching to the coalescence 
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FIG. 2. (colour online) PFC simulation of clustering phe- 
nomena on a system of 256x256 atoms after 225,000 time 
steps containing clusters with various sizes and concentra- 
tions; (a) The developed structure of a long-lived cluster; (b) 
The initially distorted structure; For graphical illustration, 
the concentration field is superimposed on the density field, 
and ranges from dark blue to dark red as the Cu-content in- 
creases. 



of excess vacancies. They have further suggested that 
the solute atoms segregate towards the loops stabilizing 
them into solute clusters. Also, tracing vaca ncy clus- 
ters by positron annihilation, Somoza et ah'sHSI have pro- 
posed that vacancy-Cu pairs are present at the quenched- 
state in Al-1.74at.%Cu alloy. To our knowledge, our 
PFC simulations are the first atomic-scale simulations 
to support the above hypothesis of vacancy/dislocation- 
mediated solute clustering and nucleation mechanisms of 
early stage precipitation. 



B. Analysis of work of formation 

We further investigated the above mechanisms of clus- 
ter formation and growth by analyzing the system en- 
ergetics for a long-lived cluster. To avoid possible finite 
size effects, a test with same conditions as those of the 
above simulation was performed on a larger system, e.g., 
512x512 atoms. The strain field caused by the disloca- 
tions displacement fields is evaluated by 

Ntr. 3 . s 

i—1 j—1 ^ ^ 
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, which is calculated over triangulated density peaks us- 
ing the Delaunay Triangulation method. Ntri is the num- 
ber of triangles in the field, Oq is the dimensionless equi- 
librium lattice parameter (the number of grid points re- 
solving one lattice spacing, i.e., 8), aij is the length of the 
j*'' side of the triangle. Small clusters, each accom- 
panied by at least one dislocation, appear to be in local 
equilibrium with the matrix shown in Fig. P[a)). Dur- 
ing the simulation, following Fig. and (c), cluster 
"a" continues to grow while, simultaneously, its accom- 
panying dislocation climbs up towards nearby disloca- 
tions, creating larger local strain fields (i.e. 0.001, 0.0016 
and 0.014 for cluster "a" in Fig. ^&), (b) and (c), re- 
spectively). This mechanism of stress relaxation through 
solute segregation has been shown through phase-field 
studies by Leonard and Haataja'^^ to be the main cause 
of alloy destabilization by structural spinodal decompo- 
sition in the presence of dislocations. Also, PFC studies 
of thin layers deposition by Muralidharan and Haataja^"^ 
indicated that, due to the above mechanism, some immis- 
cible alloys exhibit miscibility gap around the inter-layer 
interface in the presence of coherency stresses. 

The effect of dislocations on the nucleation of clusters 
is investigated by considering the following form of work 
of formation: 

H/ = 27ri?7 + 7ri?2(-A/-f AG^)- AG^^-f AGd (7) 

where R is the cluster radius in terms of number of lattice 
spacings and 



7 = 



J Area I I 



(8) 



is a Cahn-Hilliard type interfacial free energy per unit 
length of the interface in 2D. Area represents the area of 
the surface containing the cluster, and L is the circum- 
ferential length of a round cluster of radius R. Assum- 
ing low dislocation density in the system, the interfacial 
free energy is taken to be solely chemical, neglecting the 
structural contributions^. 



(9) 



is the bulk driving force for nucleation of a cluster at a 
given concentration, where superscripts '6' and 'cr de- 
note the bulk matrix and cluster "phase" quantities, re- 
spectively. 



AG, = 2GaS^ 



Ga 



(10) 



represents the strain energy for a coherent nucleu^^H^ 
where S is the misfit strain and Ga and Kb are 2D shear 
and bulk moduli, respectively, calculated from PFC 2D 
mode approximation ■ . 



AGsr^v\dEAln{R) 



(11) 



is defined as the stress relaxation term due to segrega- 
tion of solute into dislocations^, where A — -^^^z^, 



V = — 1, 77 = is the linear expansion coefBcient 
with respect to concentration, E is the 2D Young's mod- 
ulus^'^^, Xd = (f^)~^, ^bf represents a weighted average 
of the burger's vectors around the dislocations accom- 
panying the cluster and a is the lattice parameter. The 
prefactor of the logarithm term, rf'XdEA, approximates 
how strain energy is reduced due to solute segregation 
around a dislocatioil23. 



AGd = CA 



(12) 



accounts for the increase in the total system energy due 
to presence of dislocations, where ^ is a prefactor of or- 
der ten givi ng the average amount of energy per disloca- 
tion core^. Fig. [sj^d) plots the evaluation of the above 
form of work of formation (Eq. ([7|) for cluster "a" at 
different mean concentrations up to that of the largest 
cluster shown in Fig. [sjjc). The mean concentration of 
each cluster is estimated within a radius of R, defined by 
radially averaging the radius of the concentration field 
bound by a threshold of [c^ 



N 



The dashed 

curve in Fig. [3|^d) represents the work of formation for 
direct homogeneous nucleation of clusters in absence of 
dislocations, i.e., J^bf = 0. The energy barrier for ho- 
mogeneous nucleation seems to be smaller than that of 
the dislocation-assisted clustering by a single dislocation, 
i.e., E6| = 1. However, according to the plots shown in 
Fig. [sjjd), as CahiPH also pointed out, the barrier for 
formation of clusters on dislocations can be significantly 
reduced or even completely eliminated by increasing the 
magnitude of strain field around the dislocations (i.e. in- 
creasing S&f ). Notably, the local minimum also shifts to 
larger nucleus sizes until it vanishes (i.e., work of forma- 
tion continuously slops down vs. R). 

It is noteworthy that, in the absence of quenched-in de- 
fects, nucleation of the second phase requires introduc- 
tion of a thermally-activated noise to produce fluctua- 
tions in both density and concentration fields. Assuming 
dislocations are present in the bulk matrix of a supersat- 
urated quenched alloy, in this study, we demonstrate how 
elasticity itself can drive the system into a phase transi- 
tion. The infiuence of a thermally-activated noise on the 
transformation kinetics will be investigated in a future 
study through use of a well-defined noise algorithm. We 
have, however, observed in our simulations that in the 
case of a mismatch between the two species, such as in 
Al-Cu alloys, introducing a Gaussian noise to both den- 
sity and concentration fields will not have a major impact 
on the overall path of the transformation. In other words, 
the phase transformation is mainly driven by the inter- 
actions between the elastic fields of the dislocations and 
the solute atoms. 

The total work of formation, AGtot, is also estimated 
numerically by measuring the change in the grand poten- 
tial within a box engulfing cluster "a" during its forma- 
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FIG. 3. (colour online) (a-c) Snapshots taken at 3 difTerent times showing the structural changes during formation of cluster 
"a"; (d) work of formation (evaluated from Eq.( [t])) vs. R for increasing dislocation strain fields, i.e., increasing Efe?; the 
dashed curve represents the work of formation for direct homogeneous nucleation of clusters in absence of dislocations, i.e., 
when = 0; (e) the variation of numerically evaluated total energy, AGtot, and weighted average burger's vectors, Eb^ , due 
to the formation of cluster "a" in the above box; 



tion and growth in the bulk matrix, i.e., 

AGtot = f n~ f n'' = 

Jv Jv 

f [f-^i,c-^l„n]- f [f - ^ilc'' - f.in"] (13) 

Jv Jv 

. Here, ^J'c = '^^ and /i„ — ^ are diffusion potentials of 
concentration and density fields, respectively, and V is 
the total volume. The above work of formation has con- 
tributions from the interfacial energy and driving force 
for formation of clusters (i.e., AGtot = AG^—AGv), both 
of which include the elastic effects. Since the above box 
contains only one cluster, the calculated change in the 
grand potential accounts for the structural and composi- 
tional changes during the formation and growth of only 
cluster "a". While the growth of cluster "a" raises the 
local free energy, other parts of the system may undergo 
a process of annihilation and/or shrinkage of sub-critical 
clusters and their accompanying dislocations leading to 
an overall decrease in the free energy of the system. As 
can be seen in Fig. [sjje), the total work of formation in- 
creases with the growth of cluster "a" until a maximum 
value, after which it starts to decrease. Also, as can be 
seen in this figure, the estimated values of 'Sbf at vari- 
ous sizes of cluster "a" closely corresponds to its analyt- 
ical relationship with the cluster size at the local min- 
ima mapped on the energy plots of Fig. [sjjd). Likewise 



the work of formation, during formation and growth of 
cluster "a" , the value of "Ebf reaches a maximum at the 
critical size of the cluster. 

According to our data, cluster "a" continuously grows 
in presence of dislocations implying that, at each sub- 
critical cluster size, the system is sitting at a local energy 
minimum. Since cluster "a" at each sub-critical size is in 
a local equilibrium with the matrix we call it a metastable 
precursor to the cluster "a" with a critical size. This is 
analogous to previous PFC studies of crystals solidifica- 
tion which show that metastable amorphous precursors 
emerge first due to their lower nucleation barrier than 
that of a crystalline solicp2l3ni_ jj^ our case, the nucle- 
ation barrier is lowered by the effect of locally straining 
a sub-critical cluster (as a result of local accumulation of 
dislocations burger's vectors, as illustrated in Fig. [2][a)), 
making it thermodynamically favorable for the cluster to 
receive more solute atoms from the matrix and grow in 
size. 



C. Metastable phase coexistence 

The metastable coexistence between a sub-critical clus- 
ter "a" and the matrix at the quench/aging temperature 
is elucidated by evaluating the mean field free energy 
of a system comprising an unstrained matrix phase and 
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FIG. 4. (colour online) Common tangent construction using 
mean-field free energy curves of unstrained (solid curve) and 
strained solid phases (dashed curves). 

strained solid phases with different magnitudes of distor- 
tion, i.e., a uniform strain. The free energy-concentration 
curve of a strained solid phase, at a given temperature, 
can be achieved by calculating the peaks of correlation 
kernel , at locations slightly off those of the equilib- 
rium density peaks, kj, for a square structure. The in- 
troduced amount of strain is defined in Fourier space by 

e^\k~k,\/k, (14) 

, where index j denotes one family of planes in recipro- 
cal space. As can be inferred from Fig. |4j increasing the 
amount of strain from 0.0016 to 0.014 (corresponding ap- 
proximately to the average strain within the cluster "a" 
shown in Fig. [sjjb) and (c), respectively), raises the free 
energy in the strained solid. The free energy wells also 
shift to different concentrations of solute. Such a configu- 
ration admits a common tangent between the free energy 
curves of the unstrained matrix (e.g. the solid curve) and 
the distorted ones (e.g. dashed curves)P3, leading to a 
(metastable) multiphase coexistence with a lower free en- 
ergy (as demonstrated in Fig. |4]). In other words, at each 
level of local strain, there is a thermodynamic driving 
force for a transformation from a single-phase structure 
of a strained matrix to a phase-coexistence between a 
strained cluster and an unstrained matrix. On the other 
hand, despite the fact that the above transformation is 
thermodynamically favorable, the configuration of energy 
plots in Fig. [4] implies that the driving force for nucleation 
is lower for the strained cluster (i.e., using the definition 
of A/ in Eq. ([7|). However, since the energy curves in 
this figure are derived from a mean-field PFC approxima- 
tion, the illustrated phase-coexistence does not carry the 
effect of interfacial energy and only includes a mean-field 
sense of the misfit strain. These factors have a signifi- 
cant impact on the thermodynamics of phase-coexistence 
at cluster sizes smaller than that of the critical nucleus. 



In fact, the previously described stress relaxation term 
in the definition of work of formation (Eq. (|Tl|)), AGsr, 
overcompensates for the effect of reduced driving force 
for formation of strained clusters. 



D. Clustering mechanism 

Based on our PFC simulations, we propose the fol- 
lowing mechanism of clustering: (1) Stress relaxation 
by segregation of solute atoms into highly-strained areas 
in the matrix, such as around dislocations, (2) strain- 
aided nucleation of sub-critical clusters at concentrations 
higher than that of the matrix and (3) subsequent growth 
and enrichment of sub-critical clusters into overcritical 
sizes, only if a sufficient strain field is preserved, to over- 
come the nucleation barrier. The above mechanism is 
consistent with the experimentally observed formation 
and enrichment of highly-strained coherent GP-zones in 
quenched-aged dilute Al-Cu alloy^, proposed as the ini- 
tial step before precipitation of the semi-coherent and 
incoherent equilibrium 0-phase^. GP-zones in dilute 
binary Al alloys are normally known as coherent/semi- 
coherent particles often with a crystal structure and com- 
position similar to those of the final equilibrium precipi- 
tat d^ l ^^l Our clusters possess the same chemical compo- 
sition and lattice parameter as those of the equilibrium 
theta-phase pre-set by the relevant peaks in our correla- 
tion functions. Thus, they would represent an early-stage 
evolution of the so-called GP-zones. An investigation 
on the transformation of GP-zones into the subsequent 
metastable and equilibrium precipitates will be followed 
in a future study in 3D with more complex crystal struc- 
tures. We expect to observe a gradual loss of coherency 
as GP-zones grow in size, as dictated by the energy argu- 
ments. We also note that we expect our results to hold 
qualitatively in 3D, since the same type of elastic effects 
are expected to appear around the dislocations regard- 
less of their dimension and any possible partial splitting 
of dislocations around the clusters. 



V. SUMMARY 

In summary, we showed that the alloy phase field 
crystal model of ref.l^ which stabilizes different crys- 
tal structures can be used to simulate and analyze the 
mechanisms of clustering phenomenon in bulk lattice of 
quenched/aged alloys. In accordance with the existing 
experimental observations, our simulations suggests that 
quenched-in defects, such as dislocations, significantly 
lower the energy barrier for nucleation of clusters. Fur- 
thermore, analysis of overall system energy and local en- 
ergy changes reveal that the formation and growth of 
sub-critical clusters are thermodynamically favorable in 
conjunction with quenched-in mobile dislocations. Con- 
sistent with existing experiments, our simulations shed 
significant light on the elusive energetic mechanism of 
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